(* ::Package:: *)

(*Version 6
Author: Ian Brubaker
CoAuthors: Dr. Darrel Smith & Dr. Tim Callahan*)
BeginPackage["DetFit`","ErrorBarPlots`"]
DetFit::usage="DetFit[{\!\(\*SubscriptBox[\(x\), \(1\)]\),\!\(\*SubscriptBox[\(x\), \(2\)]\),...,\!\(\*SubscriptBox[\(x\), \(n\)]\)},{\!\(\*SubscriptBox[\(y\), \(1\)]\),\!\(\*SubscriptBox[\(y\), \(2\)]\),...,\!\(\*SubscriptBox[\(y\), \(n\)]\)},{\!\(\*SubscriptBox[\(\[Sigma]\), \(1\)]\),\!\(\*SubscriptBox[\(\[Sigma]\), \(2\)]\),...,\!\(\*SubscriptBox[\(\[Sigma]\), \(n\)]\)}, \!\(\*
StyleBox[\"domain1\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"domain2\",\nFontSlant->\"Italic\"]\), \!\(\*
StyleBox[\"options\",\nFontSlant->\"Italic\"]\)] Uses the determinates of equivalent-length-lists for x, y, and weighted uncertainty, \[Sigma], values to plot a linear regression of the given data. Domains 1 and 2 determine the domain of the fit line and must be real numbers. Options are the same as the options in Show, which does not includes options for plot functions."
DetFit::lengthmismatch="The arguments x, y, \[Sigma] must be of equal length."
Begin["`DetFit`"]
DetFit[x__List,y__List,\[Sigma]__List ,r1_,r2_,opts___(*underscored make the variables a function argument*)]:=
If[(Length[\[Sigma]]==Length[x])&&(Length[x]==Length[y]),Module[{delta,a1,b1,deltaa,deltab,cm,cr(*these are the "local vaiables" used only in the function and not used for an output*)},
(*The following comes directly from Dr. Smith's code, which can be referenced on his website physicsx.pr.erau.edu*)
delta =Sum[(1/\[Sigma][[i]]^2),{i,Length[x]}]* Sum[(x[[i]]^2/\[Sigma][[i]]^2) ,{i,Length[x]}]- (Sum[(x[[i]]/\[Sigma][[i]]^2),{i,Length[x]}])^2;   
a1 = (1/delta)*(Sum[(x[[i]]^2/\[Sigma][[i]]^2),{i,Length[x]}] *Sum[(y[[i]]/\[Sigma][[i]]^2),{i,Length[x]}] - Sum[(x[[i]]/\[Sigma][[i]]^2),{i,Length[x]}]* Sum[(x[[i]] y[[i]]/\[Sigma][[i]]^2),{i,Length[x]}]);
b1 = (1/delta)*(Sum[(1/\[Sigma][[i]]^2),{i,Length[x]}]* Sum[(x[[i]] y[[i]]/\[Sigma][[i]]^2) ,{i,Length[x]}]- Sum[(x[[i]]/\[Sigma][[i]]^2),{i,Length[x]}]* Sum[(y[[i]]/\[Sigma][[i]]^2),{i,Length[x]}]);
deltaa = Sqrt[(1/delta)*Sum[(x[[i]]^2/\[Sigma][[i]]^2),{i,Length[x]}]];
deltab = Sqrt[(1/delta)*Sum[(1/\[Sigma][[i]]^2),{i,Length[x]}]];
cm=Sum[((y[[i]]-(b1*x[[i]]+a1))/\[Sigma][[i]])^2,{i,Length[x]}];
cr=cm/(Length[x]-2);
Column[{Manipulate[Show[ErrorListPlot[Table[{{x[[i]],y[[i]]},ErrorBar[\[Sigma][[i]]]},{i,Length[x]}],PlotStyle->PointSize[size]],
Plot[b1*f+a1,{f,r1,r2}],opts],{{size,.012,"Point Size"},.01,.15}],Print["Intercept  a =  ",a1, "  \[PlusMinus] ",deltaa] ,Print["Slope      b = ",b1, " \[PlusMinus] ",deltab],Print["\!\(\*SubscriptBox[SuperscriptBox[\(\[Chi]\), \(2\)], \(Minimum\)]\) = ",cm],Print["\!\(\*SubscriptBox[SuperscriptBox[\(\[Chi]\), \(2\)], \(Reduced\)]\) = ",cr]}]
(*Notice the lack of a semi-colon, this is becuase the value here will be outputted when calling this certain function*)],
Message[DetFit::lengthmismatch,"Argument \[Sigma] calls for a list of values with length equal to the number of data points"]];
End[]
EndPackage[]
